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Abstract 

An asymptotic interface equation for directional solidification near the absolute stabiliy limit is 
extended by a nonlocal term describing a shear flow parallel to the interface. In the long-wave 
limit considered, the flow acts destabilizing on a planar interface. Moreover, linear stability analy- 
sis suggests that the morphology diagram is modified by the flow near onset of the MuUins-Sekerka 
instability. Via numerical analysis, the bifurcation structure of the system is shown to change. 
Besides the known hexagonal cells, structures consisting of stripes arise. Due to its symmetry- 
breaking properties, the flow term induces a lateral drift of the whole pattern, once the instability 
has become active. The drift velocity is measured numerically and described analytically in the 
framework of a linear analysis. At large flow strength, the linear description breaks down, which 
is accompanied by a transition to flow-dominated morphologies, described in a companion paper. 
Small and intermediate flows lead to increased order in the lattice structure of the pattern, fa- 
cilitating the elimination of defects. Locally oscillating structures appear closer to the instability 
threshold with flow than without. 



PACS numbers: A7.5A.+I, 05.70.Np, 81.30.Fb, 64.70.Md 



I. INTRODUCTION 



Directional solidification is an important experimental procedure both from the applied and the 
fundamental points of view. 

On the one hand, growth techniques such as zone melting and the Bridgman method were 
developed in the fifties to purify semiconductor materials. Subsequently, the basic process became 
an important industrial tool, and it continues to remain the basis of fundamental metallurgical 
techniques. 

On the other hand, when an alloy is solidified by pushing its melt at a constant velocity along an 
imposed thermal gradient towards lower temperatures - this realization of directional solidification 
is particularly amenable to theoretical analysis - the interface between the liquid and the solid, 
which is flat at zero velocity, undergoes a morphological instability. Fundamental interest in this so- 
called MuUins-Sekerka (MS) instability motivated many current research studies of directional 
solidification. The instability is driven by impurity diffusion, it appears at a critical speed, and 
once it has set in, new structures develop, depending on a number of factors such as the interface 
roughness, the segregation coefficient, impurity diffusion, and the impurity concentration in the 
liquid. If the latter is high enough, two-phase eutectic structures may arise [g[ , which are not created 
by the MS instability. In this article, we will rather focus on the case of smaller concentrations 
leading to the growth of a dilute alloy forming a single-phase solid. 

Experimentally, pattern formation in directional solidification has been studied in great detail 
using transparent model substances |^ in order to gain fundamental insight_P,|[^,0,||]. It has also 



12| . Recently, successful 
loys have opened 



been investigated in metals, with at least a view to applications P,pl],p^, 
observations of three-dimensional directional solidification of transparent a 
the road to an understanding that goes beyond the description of two-dimensional structures, for 
which a large amount of theoretical work does exist ||l|,|l|,[|,^^|20[|2^,|2|,|2|,|2|,|2|,||,|^,||. 

As a consequence, it has now become an important task to follow up with theoretical work on 
three-dimensional directional growth. Detailed analysis of three-dimensional patterns has so far 
largely been restricted to free growth of small structures such as single dendrites, both analyti- 
cally [^,^,^ and numerically |32,33, Mj. Three-dimensional directional solidification has been 
considered in the context of phase-field modeling for small systems consisting of sonic ten cells |]35| . 
Only fairly recently, larger systems containing several thousand cells have been treated within the 
framework of an asymptotic interface equation Js^ . The main advantage of this equation is that 
it reduces the dynamical problem to a description of the interface alone, which replaces the 3D 
problem with a 2D one, and that, moreover, it is a local description, thus rendering large systems 
tractable. It was obtained by a Sivashinsky-type expansion about the point of absolute stabil- 
ity 1^^, neglecting both solute trapping and deviations from local interface equilibrium. As these 
approximations become doubtful for large solidification speeds, this model equation had to be con- 
sidered a qualitative description only and comparison with experiments had to rely on genericity 
arguments. However, we believe that the equation would be quantitatively appropriate for the 
description of directional ordering experiments on liquid crystals |^,^,|3^] and hence we suggest 
the more intriguing dynamic long-time effects predicted by it ||3^ to be observable in liquid crys 



tals, which also have response times more accessible to experiments 40 1. In any case, statistical 
properties obtained from the simulation compared reasonably well with the same properties mea- 
sured in experiments. This holds apart from a systematic deviation allowing us to conclude that 
experimental cellular patterns do not start growing from a random (Poissonian) distribution. 
In earth-bound experiments, convection always exerts a major influence on the arising patterns 
Therefore, we extend our previous work here to include a flow in the description. Since this is 
but a flrst approach to a large-scale system with flow, in order to facilitate the subsequent analysis 
of results, we consider the conceptually simplest case, an externally imposed flow parallel to the 
interface, approaching a constant velocity far from it. For this situation, representing the "poor 
man's convection" ||4^, an asymptotic description has been derived by Hobbs and Metzener p3[ |. 
One-dimensional interfaces in the presence of flow have been studied in ji^ . Whereas in the latter 
investigation solute trapping and the deviations from equilibrium of the interface were fully taken 
into into account using the Aziz model [ ^St , we will restrict ourselves here to the considerably 
simpler situation of a constant segregation coefflcient and neglect interface kinetics, to separate 
generic effects from those that are only important at high speed. Also the full equation is more 
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susceptible to numerical instabilities, because its solutions tend to develop large amplitudes taking 
them out of the validity domain of the expansion The only difference between the interface 
equation employed in |pq] and here will therefore be the term generated by the flow. It renders 
the equation nonlocal which already makes its simulation much more tedious than that of the 
original equation and restricts our simulations to a few hundred solidification cells rather than a 
few thousand. 

This article consists of five sections. In Sec. II, we give the basic equations together with a linear 
stability analysis, leading to the prediction of new morphologies induced by the flow. Section III 
describes the numerical approach and our methods of velocity measurement. We give a number of 
simulation results in Sec. IV, exhibiting the basic patterns in the part of parameter space where 
some degree of spatial order prevails. A more detailed analysis and characterization of these struc- 
tures is presented in the subsequent article |^6|| , henceforth referred to as II. Conclusions regarding 
flow effects on pattern dynamics and the bifurcation structure of the system are summarized in 
Sec. V. An appendix provides the real-space representation of the flow term thus displaying its 
nonlocal nature. Some analytic approximations are worked out in detail in a second appendix. 



II. MODEL EQUATION AND LINEAR STABILITY ANALYSIS 



A. Long-wave equation 



The equation to be studied is 



Ctt - (2 + i + \/\t + (1 + ^ + v'c + sfcv^c + sfcGc - f-^m 

= 2CtS/\ + 2 (|VC|')^ - (|VC|') - 2:.(VC)V (v^c) - {(vov^c} 

-2V{(VC)|VC|'} . (1) 

We denote by V the two-dimensional gradient operator, V = {d^, dy). The position of the liquid- 
solid interface is ^(x, y, t) and coordinates are chosen such that the temperature gradient is oriented 
along the z axis, whereas the flow is parallel to the x axis. C[C] is a linear but nonlocal functional 
of given via its Fourier transform 

^[£[C]](p)-|p|.F[C](p), (2) 

with !F denoting the transform, — J^^dx J^^dy ({x) exp{—ipx), and p its wave vector 

argument. Position space representations for C[Q in one and two dimensions are given in App. A. 

There are four nondimensional parameters in our equation: the segregation coefficient k, the ratio 
p = Dg/D of the diffusion coefficients for impurities in the solid and the liquid, the nondimensional 
temperature gradient G, and the strength / of the flow. The one-sided model is characterized by 
1^ = 0, the symmetric model by 1/ = 1. G is related to physical parameters via: 

- ^ SD'L^mAc G 

In this expression, L is the latent heat (per unit volume) of the transition, m the absolute value 
of the slope of the liquidus line, Ac the miscibility gap at the base temperature Tq (Tq = — 
TOCo, with Tm the melting temperature of the pure substance, cq = Coc + Ac, and Coo the initial 
concentration of the liquid), 7 is the surface energy, assumed isotropic here, Va is the velocity 
corresponding to the absolute stability limit, given by Va — mLAcD/jT„ik. G is the temperature 
gradient and V the pulling (or pushing) velocity. Finally, / is given in terms of dimensional 
quantities as 

f - ^ (4) 
^ eVS ' ^ ' 
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Voo being the speed of the imposed flow far from the interface, S = Vk/ D is the Schmidt number 
{vk is the kinematic viscosity of the hquid), and e is the nondimensional distance from absolute 
stabiUty, e — l/2k — "/T„iV/2mLAcD. f is dynamically determined by the ratio of the flow and 
pulling speeds. 

Equation (|l|) is obtained from an expansion about the absolute stability limit, where the wave- 
length of the most unstable mode diverges as l/-y/e. It takes the form of a strongly nonlinear 
long-wave equation, in which wavelengths have been rescaled by y/e to make them 0(1). Thus, 
nondimensional lengths arc measured in units of a rescaled diffusion length along the x and y 
directions, with a scaling factor 1/%/e; lengths in the z direction are measured as multiples of the 
unsealed diffusion length. The time unit is a diffusion time {2D/V'^) scaled by 1/e. For fc = 0, the 
equation becomes indefinite, because there is no absolute stability limit in this case. 

For simplicity, we shall set fc = 1 in the following, i.e., Ac = const, independently of the reference 
temperature Tq. From previous experience, the choice of k is not expected to have a strong influence 
on results, as long as k does not become very small. 

To further reduce the parameter space to be explored, we also set = 1, thus restricting ourselves 
to the symmetric model. This is not a particularly realistic model for directional solidification, 
where 1^ = would be more appropriate. However, for directional ordering in liquid crystals, it is a 
better approximation than the one-sided model. Moreover, the choice v = 1 is least problematic in 
numerical simulations with a direct finite-difference discretization, because deep grooves that may 
trigger numerical instabilities are less likely to evolve with diffusion in the solid allowed. Since the 
most efficient way to deal with the nonlocal term is to work in Fourier space, we have developed a 
pseudospectral code, which is much less susceptible to these problems. Nevertheless, a comparison 
with previous results is easier, if we keep at the value then used, and experience suggests |Q 
that generic patterns are not strongly influenced by this choice. Finally, because our results must 
not be expected to be quantitative except for liquid crystal systems, we stick to the value v — 1 
here. Tests with various nonzero values of v have been performed but have not revealed interesting 
differences. 

Therefore, the important parameters to be varied in the simulations are G and /. 

The form of Eq. (|) can be guessed from symmetry and scaling arguments. The linear terms 
on the left hand side are determined (including their coefficients) by the linear stability analysis 
of the full three-dimensional model, involving diffusive transport and the coupling to the Navier- 
Stokes equations. Scaling arguments tell us that the nonlinear terms can contain at most four 
spatial derivatives and for each temporal derivative present there must be two spatial derivatives 
less (since wavenumbers scale as y/e but frequencies as e). In the absence of a thermal gradient, we 
have translational symmetry in the z direction, so we know that all nonlinear terms must contain 
derivatives of ^ only. If there is no flow, we also have parity symmetry, which constrains the number 
of spatial derivatives of terms not containing / to being even. Finally, rotational symmetry can 
be invoked to exclude terms such as (V^)^ ||2^]. From these considerations, one obtains all the 
nonlinear terms on the r.h.s., but not their prefactors, of course, for which the full expansion must 
be performed What can also be guessed is that the flow should lead to a nonlocal term 

breaking the parity symmetry with respect to the x coordinate. It is not clear beforehand, however, 
that it does not introduce additional nonlinear local terms. The nonlocal term on the left-hand 
side is, in a sense, the simplest nonlocality possible (see Appendix A). 

B. Linear stability analysis 

The problem of coupled morphological and convective instabilities has a long history of detailed 
study |^,^|5^,^,^2|| . In the case of buoyancy-driven convection with a lighter solute and solidifi- 
cation proceeding upward in the gravitational field (leading to unstable stratification), it was found 
that in most cases the coupling between the instabilities is weak due to a large disparity in unstable 
wavelengths |5^ , even though an oscillatory instability may occur in special circumstances. 
For Rayleigh numbers below the critical value, convection delays the MuUins-Sekerka instability in 
the limit of small segregation coefficients, where a long- wave equation different from ours can be 
derived. Forced flows were studied by Coriell et al. and by Forth and Wheeler Q , who found 
that for two-dimensional disturbances the flow delays the morphological instability, whereas for 
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disturbances with wave vectors perpendicular to the flow the latter normally does not affect the 
critical conditions of the MS instability. However, for small-wavenumber modes, the MS instability 
can be enhanced |Q. Near absolute stability, the wavelength is much larger than the diffusion 
length and parallel flow was shown to be destabilizing to disturbances that travel against it | ]55[ |. 
This is the situation encountered here. 

Equation (|l|) has the steady-state solution ^ = 0. Obviously, the linear stability analysis of this 
solution will involve only the terms on the left-hand side of the equation. Inserting the perturbation 
ansatz [x = {x, y)] 



C = Ci exp {ujt + i qx) 
into Eq. (|]), we obtain the dispersion relation {q = |q|) 



OJ 



8kG -ifqxq = Q . 



(5) 



(6) 



The terms stemming from the partial derivatives are obtained in a straightforward manner, only 
the one arising from the nonlocal term may require some explanation. To compute the nonlocal 
term for the perturbation (|^) we first take the spatial Fourier transform of C, which is simply 
47r^Ci ex-p{u!t)S (q + p), then multiply it by |p| to obtain the transform of £[C]. Transforming back 
we get Ci|q|exp(wt + iqx), the derivative of which with respect to x produces a prefactor iq^- 
After dropping the common exponential factor and the prefactor of all terms, we are left with 
the q dependent expression of Eq. (^) . 

Setting Lu = LOr + iuji and decomposing the dispersion relation into its real and imaginary parts, 
we obtain 



cj^ - + ( 2 + i + 1/ ) ujrq'^ + 



The unstable mode takes the form 



1 



1 



2 I 4 

V \ q 



8kq^ + 8fcG = , 



2ljJr 



1 



2 + - + H'Z^ 



= fqxq ■ 



(7) 
(8) 



C = Cl 6Xp < ijJrt + i 



qx\xA 1] + qyy 



(9) 



i.e., it corresponds to a traveling wave moving along the x axis at velocity Vd = —'^ijqx- On the 
neutral surface, = 0, which implies 



hx 



(2 + l/A: + ^)g ' 



(10) 



hence we have a Hopf bifurcation whenever the flow is different from zero and the pattern is 
oriented such that qx ^i). The velocity of the corresponding traveling wave is simply 



/ 



{2 + \lk^v)q ■ 

Inserting ( p^ ) into (|^), we arrive at the equation for the neutral surface: 



Aq)^ 



1 

^^k 



8kG 



fqx 



{2 + l/k + i^)q 



0, 



(11) 



(12) 



where the last term depends on the angle between the flow and the wave vector only, not on the 
modulus of the latter, a situation that we describe by setting cos0 — qx/q- Obviously, for fixed 
values of / and (/) the essential change of the neutral curve (in the qG plane) brought about by 
the fiow is an increase of the critical value of G where the instability first appears. Hence the 
flow has a destabilizing effect, since the region of parameter space, where the planar solution is 
unstable corresponds to values of G below the critical value. For given G, the flow can be made 
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large enough to render the planar front unstable even near q = 0, i.e., with respect to homogeneous 
perturbations. Figure |l| displays the function g(q) for two values of the flow (assuming cos(/) — 1). 
With weak flow, g{q) is positive at q — 0, and if G is small enough, the function has two zeros at 
positive g, i.e., a band of modes containing only finite q values is unstable. For strong flow, g{q) is 
negative at g = 0, i.e., all modes up to the marginal one are unstable. 




0.0 0.5 1.0 1.5 2.0 

q 



FIG. 1. The function g{q) determining the neutral surface. G = 0.5. Upper curve: f — 0, lower curve: 
/ = 10. 

The neutral surface is the set of zeros of g{q) in the space spanned by g, G, and /. For 
convenience, we restrict ourselves to the qG plane and draw neutral curves for several discrete 
values of the flow parameter / in Fig. ||, to demonstrate both the G and / dependences of the 
neutral surface. 




0.0 0.5 1.0 1.5 

q 

FIG. 2. Neutral curve for different flow strengths /. Lower solid line: f — 0, dashed line: / = 1, 
dash-dotted line: / = 2, dotted line: f ~ 3, solid line: / = 4, higher dashed line: / = 5. The curves shift 
upward with increasing / by an amount proportional to 

A calculation of the critical value of the temperature gradient in the general case (requiring 
duJr/dq = 0) leads to 
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1 

1 + k + 8{l + 2k + ki^y 



exhibiting the fact that the instabihty threshold depends on the angle between the flow and the wave 
vector of the perturbation. Without flow, the threshold is given by the first term of Eq. (|l3|) |Q 
and the bifurcation, transcritical in two dimensions, is known to give rise to hexagonal patterns 
at onset fs^ . From the equation, we can immediately conclude that for values of the temperature 
gradient satisfying G > Gdf = 0), there exists a critical flow strength given by 



2 _ 8{1 + 2k + kiyf ( - 2fc2 



above which the planar front is destabilized by the flow alone. In this case, the patterns emerging as 
the planar interface becomes unstable will not be hexagons but rather stripes oriented orthogonally 
to the flow (since these are the most unstable disturbances). They should drift against the flow 
with a speed approximately given by Eq. ( pi] ) and calculated more precisely below. As we shall see 
in Sec. IV, these predictions are borne out by the simulations. It is then an interesting question, 
how the system will behave on decrease of G below the zero-flow threshold. What happens when 
the pattern amplitude approaches saturation cannot be predicted from the linear analysis but will 
be discussed in some detail in II. 

To get an idea of the behavior of the drift speeds to be expected beyond the bifurcation point, 
we compute the fastest-growing unstable mode, which should provide a decent approximation 
to observed wavelengths close to threshold. For simplicity, this calculation is restricted to the 
symmetric model {v — 1) with k ~ \. The q value, at which the growth rate is maximum is 
obtained by differentiating Eqs. (^, (|^) and setting dujr/dq ~ 0, which gives us two more relations 

-2w,^ + 8(7Wr + 12g^-16g = 0, (15) 
dq 

8quj, + ^ {2ur + V) = , (16) 

from which the four unknowns w^, diOi/dq and q can be determined. Two simplifications are 
straightforward, giving expressions for uji and dcoi/dq: 



dq 2ujr + 4g2 
flxq 

= ^ 

2uJr + 4g2 



{2fq^ ~ Sqto,) , (17) 

(18) 



which leaves us with two equations for ujr and q. 

Incidentally, we can immediately gather an interesting consequence from these equations regard- 
ing the question of convective versus absolute instability p7| . If we require w,- = in Eq. (p^, 
this implies duji/dq = by virtue of Eq. ([1^). Hence, at the critical point of the linear instability, 
we have dio/dq = 0, i.e., the group velocity of a localized perturbation vanishes. This means the 
thresholds for convective and absolute instabilities coincide in our system, which therefore is never 
only convectively unstable. This statement remains true for arbitrary values of k and v. 

In the following, we will assume that the stripe pattern is oriented orthogonally to the flow (as 
it usually is if it arises spontaneously from a random initial condition), therefore qx = q. The 
equations determining the two remaining unknowns are then 

£2 4 

'^r - + 4t^r9^ + Sg" - + 8G = , (19) 



Pq'^LUr 



2{uJr + 2(72)3 



4wr - eg-" + 8 = . (20) 



The limiting cases / ^ 1 and / ^ 1 of this system of equations can be treated analytically, 
detailed expressions are given in App. B. For small flow velocities the interface drifts at a speed 
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that is proportional to /, whereas for large velocities, the drift speed becomes proportional to \fj . 
The imaginary part of the growth rate u)i is proportional to / in both cases but with different 
proportionality constants. A numerical solution of the system (|l9|), ( |20| ) is not difficult. The 
resulting "drift frequencies" a;^ are given in Fig. ^ for several pertinent temperature gradients, 
together with the asymptotic analytic expressions from App. B. We will compare theoretical with 
measured drift velocities in Sec. IV. 




FIG. 3. Imaginary part LOi of the complex growth rate lo for the fastest-growing mode ("drift frequency") 
as a function of the flow parameter / for G = 0.1 (solid line), G = 0.35 (dashed line), and G — 0.6 
(dash-dotted line). The thin dotted line is the asymptotic approximation for large /, Eq. (B19), whereas 
the other thin lines denote asymptotic expressions for small / and different G values, Eq. (|B8[). 



III. NUMERICAL APPROACH 



A. Discretization 

Equation (|l|) was simulated with periodic boundary conditions on quadratic grids of sizes between 
32x32 and 512x512. The mesh size h was usually 0.5, although at thermal gradients below 0.35, 
we had to reduce it to keep the code stable. With a 128x128 lattice, hexagonal structures would 
contain on the order of 100 cells for h = 0.5. Previous simulations have shown that this is 
roughly the size of a typical dynamical grain, inside which a system at not too low a temperature 
gradient manages to get rid of all its topological defects and to attain complete hexagonal order. 
In order to have several grains in the numerical box, a number of 256x256 systems were simulated. 

Temporal discretization was done by a simple explicit first-order Euler scheme. Two variants 
of the code were implemented; in the first, spatial derivatives were approximated by second-order 
accurate symmetrical stencils, whereas in the second, a pseudospectral approach, derivatives were 
computed via fast Fourier transform, i.e., they were accurate to order for a mesh size h and 
linear grid dimension N. The flow term was always evaluated via its Fourier representation, as 
a real-space calculation would have required the computation of a double integral on the whole 
system at each lattice point (see App. A). We will refer to the first, less accurate approach as the 
mixed code and to the second as the (pseudo)spectral one. 

The mixed code was useful for the treatment of larger systems (256x256) over longer times, since 
it was faster than the spectral code by a factor of about six in this case. However, the simulation 
of systems with smaller temperature gradients {G < 0.35) necessitated the use of the spectral 
code, both for accuracy and stability reasons. In our previous study of three-dimensional rapid 
solification ||3^ , gradients below 0.35 remained essentially inaccessible due to numerical instabilities 
arising at a mesh size of 0.5. Reduction of the mesh size mitigated the problem, but restricted 
accessible system sizes. 
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B. Velocity measurement 



In the presence of flow, patterns move laterally, so it became desirable to measure tlieir velocity. 
This was done via a correlation function method as follows. The dynamic evolution of the interface 
was simulated over a time interval extending from t to t + At. Then the quantity 

c{Ax,Ay,At,t) = (^[C{x + Ax,y + Ay,t + At)-ax,y,t)f) (21) 

was evaluated for a number of values Aa; and Ay that were small multiples of the grid spacing 
h. Angular brackets denote spatial averaging (over the whole grid). Next the minimum of the 
correlation function c was determined via parabolic interpolation from the surroundings of the 
minimum value obtained within the discrete set {Ax, Ay}. An approximation to the velocity at 
time t + At was then obtained as v = {vx,Vy) with Vx — Ax* / At and Vy = Ay*/ At, where Ax* 
and Ay* were the coordinates of the minimum. 

We tested this procedure on a variety of analytically prescribed interfaces. It turned out highly 
reliable and accurate (in the ppm range and better) whenever the interface had a constant shape, 
its only dynamics being a lateral drift motion, and the time step At was not chosen too short. 
The accuracy deteriorated to fall into the percent range, when shape changes were allowed. This is 
understandable, as with a shape-changing interface the drift velocity is not even precisely defined. 
A one-dimensional example will clarify this point. Consider 

({x,t) = sin k{x — vt) cos ujt , (22) 

where intuitively one would associate the velocity v with the motion of the pattern. But we also 
have 

C(x,t) = isinfc X — f w — ^) t -fisinfc x—(v+'^\t , (23) 
2 L V k / i 2 L V fc/J 

that is, the pattern is decomposable into two waves drifting at different velocities v — uj/k and 
V + Lu/k. 

In this case, the correlation function can be calculated analytically: 

X — V At) cos uit cos uj(t + At) , (24) 

the spatial minimum of which is, for cosLot cosw(t -I- At) > 0, given by fc(Ax — vAt) — 2m: 
(n = 0,±1,±2, . . .). The example shows that for regular structures, Ax (as well as Ay) has to 
be kept smaller than a wavelength in order to avoid solutions with n ^ 0. Assuming n = 0, we 
find Ax/ At = v. This means that one obtains the intuitively expected result v for the velocity, 
whenever coswf does not change sign in the interval t € [t,t + At] (in these - rare - instances, the 
algorithm will yield Ax/ At = v ± n/kAt, i.e., the "velocity signal" will display peaks at (twice) 
the frequency cu). We conclude that in general the algorithm is reliable and robust. 



IV. SIMULATION RESULTS 



A. Some basic patterns 



The presence of a flow term considerably increases the richness of the system with regard to 
pattern formation. Whereas stable structures of the system without flow can be described in a 
summarizing fashion as more or less ordered hexagonal arrays of cells, which may be steady state 
(with some movement in the grain boundaries) or oscillatory (with phase shifts of « 27r/3 within a 
triangle of neighboring cells) or (weakly) turbulent js^, the system with flow has many more ways 
of organizing itself. Figures ^ through ^ may serve to give a first impression. Each of these figures 
displays a typical structure for a given temperature gradient at a moderate flow. 

In Fig. ^, the value of the temperature gradient is G = 0.7, i.e., larger than Gc(/ = 0), hence 
the planar interface is destabilized by the flow only. Therefore, no hexagonall cell structure can 
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develop initially (as long as the pattern is describable by the linear theory). We obtain a stripe 
structure containing defects, some of which disappear pretty fast, whereas the last few persist for 
a long time. The final evolution of this pattern up to several thousand diffusion times will be 
discussed in II. 

Figure ^ is at G = 0.6, where the MuUins-Sekerka instability is already present. It can be 
clearly seen that the flow has an organizing influence on the structure consisting of hexagonal cells: 
the dynamical grain boundaries separating differently oriented hexagonal domains try to orient 
themselves perpendicular to the flow, so they become parallel. We did not observe similar ordering 
of grain boundaries in simulations without flow, in which grain boundaries rather tend to form 
ringlike structures js^ . 




FIG. 4. Pattern at G = 0.7, / = 4.0, system size 128.0x128.0 (grid spacing h = 0.5), i.e., the lattice 
is 256x256. t = 150.0. Lengths are given in units of the (rescaled) diffusion length, times in units of the 
(rescaled) diffusion time. 




FIG. 5. Pattern at G = 0.6, / = 2.0, system size 128.0x128.0 (grid spacing h = 0.5). t = 2000.0. 
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At smaller temperature gradients, flows of moderate size do not appear to strongly perturb 
the basic structure imposed by the MS instability. The dynamics is of course different, since the 
entire pattern drifts against the direction of the flow. Comparisons with simulations without flow 
show moreover that structures display more order (after the same time of dynamical evolution and 
starting from the same random initial conditions) with flow than they do without. As the flow is 
increased, new structures develop that we discuss below. 

At G = 0.35, we find the flow to promote oscillatory structures. Figure ^ shows an example of a 
topologically ordered array of hexagons. Each cell has exactly six neighbors. The brightnesses of the 
cells indicate their different heights; white is high, black is low. Different sizes of the cells are due 
to their being in different phases of their basic oscillation. There is a phase shift of approximately 
27r/3 between neighboring cells. Phase coherence is not preserved throughout the entire array of 
cells as may be noted by comparison of the lower left and upper right parts of the pictures (not 
all big bright cells are exactly the same, and similar statements hold for the small bright cells and 
the dark cells). By animated visualization of a series of pictures, the oscillations are easily verified. 
Large cells are seen to become smaller and larger again, periodically. Small cells behave the same 
way, the only difference being a phase shift. Unfortunately, a movie cannot be transmitted in this 
media. 




FIG. 6. Pattern at G = 0.35, / = 2.0, system size 128.0x128.0 (grid spacing h = 0.5). t = 2000.0. 

From earlier work |]3^ , these oscillations are known to occur even without flow; an analytical 
discussion, not considering stability issues was given in p6[ . For large grid spacings {h — 0.5) 
we saw this dynamical state already at G = 0.4 with our finite-difference code. It was however 
observed to appear later, i.e., at smaller G, when the mesh size was reduced, and we estimated the 
bifurcation to 27r/3 oscillations to happen between G = 0.4 and G = 0.35 As it turns out, the 
spectral code with its higher accuracy does not yet produce the oscillations at G = 0.35, if flow is 
absent, but it does so in the presence of even small flows (/ = 1.0). In addition, the cellular lattice 
becomes more ordered, topological defects are eliminated more efficiently under flow. 

To finish this introductory tour through parameter space, we consider a temperature gradient of 
G = 0.25 (Fig. 1^), which was impossible to do with the finite-difference code because of numerical 
instabilities. Even with the spectral code we have to reduce the mesh size to obtain numerically 
stable results for G distinctly smaller than 0.35. 

The example suggests that patterns are generically weakly turbulent at such a small gradient, 
i.e., they show time-dependent nonrelaxational behavior. Nevertheless, some ordering influence of 
the flow may be noted even here (cells tend to align along the direction perpendicular to the flow) 
and becomes much more conspicuous as the flow is increased, to the extent of rendering structures 
regular again at larger flows. We will return to this question in II. 
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FIG. 7. Pattern at G = 0.25, / = 2.0, system size 76.8x76.8 (grid spacing h = 0.3). t = 1500.0. 

B. Measured properties 

In order to characterize the bifurcation from the planar front, we introduce the standard deviation 
A = y/ ((C — (C))^) as a measure for the amplitude of the steady state structure, reached in the 
late stage evolution of an initially sinusoidal interface. Figure ^ shows the amplitude so obtained 
as a function of the temperature gradient. 
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FIG. 8. Amplitude A of steady-state stripe patterns for a fixed flow of / = 2.0 as a function of the 
temperature gradient. The initial condition was a sinusoidal pattern with wavenumber 0.98, with the wave 
crests oriented parallel to the y axis. Measurement of the amplitude was done after t — 1000.0 in units 
of the (rescaled) diff'usion time. The dotted line gives the theoretical position of the bifurcation point, 
according to Eq. (^3|). 

Due to the small system size (16.0x16.0), it was possible to keep stripe structures stable well 
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below the threshold of the appearance of hexagons ((5 = 2/3). As expected the bifurcation is 
supercritical. 

In the presence of any nonzero flow, the basic structures appearing at the instability threshold 
are stripes rather than hexagons. Whereas our large-flow predictions, discussed below (and, in 
more detail, in II), might require some effort to be realized experimentally, this result should be of 
immediate experimental relevance, since it is valid no matter how small the flow. Extremely small 
flows would make the range of temperature gradients in which stripes dominate over hexagons very 
small, but this range would nevertheless be present and necessarily observed, at least temporarily, 
on crossing the bifurcation threshold via reduction of G. Since the argument for the prevalence of 
stripes is drawn from linear stability analysis, it need not continue to hold, once amplitudes become 
large enough for nonlinearities to play an important role. In II, we will see that this indeed happens 
as long as the flow is not too strong. The transcritical nature of the bifurcation to hexagons will 
then turn out important. 

Measuring the velocity of the interface for values of G which are sufficiently close to the instability 
threshold we find that it exhibits damped oscillations. Examples for G = 0.7 and G = 0.6 are 
presented in Fig. ||. Our discussion of the velocity measuring procedure in Sec. Ill suggests that this 
phenomenon may be due to some dynamics superimposed on the drift motion. This is corroborated 
by examining the frequency of oscillation. 
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FIG. 9. Oscillations of the lateral pattern velocity after initialization with a random structure, 
line: G — 0.6, dashed line: G = 0.7. For the smaller value of G, the oscillations decay faster. 



Solid 



The figure shows clearly that the frequency is slightly higher for the larger value of G (for 
G = 0.7, there are 17 oscillations in the time window displayed, but only 16 for G — 0.6). A 
precise determination of the (angular) frequency reveals that it is very close to the (angular) 
frequency uj = VSkG of homogeneous solutions |Q to Eq. (^. As has been noted before psj , 
patterns initiated close to the threshold Gc oscillate as a whole before settling down into a steady 
state. 

The oscillatory velocity pattern is thus an effect of the temporal modulation of the pattern 
[similar to the cosojt term in Eq. (|2^)] and we should consider the average over these oscillations 
the true velocity. They are a nongeneric feature of the current amplitude equation (|^), which is 
not shared by other equations such as the Kuramoto-Sivashinsky equation. For lower values of G, 
these oscillations are not present. 

Let us now look at the drift velocities measured for different values of the temperature gradient 
and the flow. Figure |l^ collects some data corresponding to gradients between G = 0.35 and 
G ~ 0.7. For flow strengths below |/| — 6, all the data points collapse approximately onto 
one curve. This is to be expected from the linear-stability result (B£), which shows that the G 
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dependence of the drift velocity is weak for small /. The dashed line in the figure shows this result 
for G = 0.6. 

On the other hand, the spreading of the data points for / > 6 does not follow from the cor- 
responding result (B20) for large /. Clearly, the linear theory breaks down here. As we shall 
see below, the deviation of the drift velocity from the analytic result happens in the vicinity of 
the transition to a new pattern, where the drift velocity is determined by nonlinear effects. We 
will discuss this in particular for one temperature gradient, but there seems to be a morphology 
transition in all cases where a strong deviation from the linear theory arises; the transition is not 
always to the same new pattern however. (For G = 0.4, there seem to be two transitions.) 
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FIG. 10. Measured drift velocity as a function of the flow for different temperature gradients. Stars: 
G = 0.7, squares: G = 0.6, triangles: G = 0.5 (most of them covered by other symbols), inverted triangles: 
G = 0.4, circles: G = 0.35. Dashed line: Analytic velocity from linear stability analysis for G = 0.6. 



Without flow, the equation of motion (|l|) is symmetric under the parity transformation x —x. 
Drifting patterns arise, because this symmetry is broken by the flow term, meaning that we do 
not have spontaneous symmetry breaking as, e.g., with parity breaking patterns in the purely 
diffusive case |56 . Whereas those patterns are not stable in extended systems, the present drifting 



structures are robust. 

It is easy to see that flow-induced drifting cells must themselves be asymmetric with respect to 
a mirror plane parallel to the yz plane, even though this asymmetry may be barely perceptible to 
the eye. To show this, we assume the opposite to be true, i.e., C to be symmetric under x ~x. 
Transforming Eq. (|^) to a comoving frame, which results in Ct ^ Ct — vCx, we obtain from the 
antisymmetric part of the equation: 



2Mivcr 



(25) 



But this relation must be invariant under an exchange of x and ~x and a replacement of v by —v. 
Hence, the flow term must vanish, which means it cannot play any role. This is in contradiction 
with the fact that drifting patterns on large scales are not observed in the absence of flow. 

We have probed the asymmetry of cells in a number of different ways as a test of our numerical 
procedure. The most distinctive method was as follows: first the 2D interface was cut by a straight 
line parallel to the x axis through the middle of a cell, to get its profile. Then the top part of the 
cell (everything above a chosen threshold height) was fitted to a parabola. The difference between 
the fit function and the actual cell profile is given for one case without and one with flow in Fig. 
Evidently, this difference is symmetric in the former and asymmetric in the latter case. Since the 
flow comes from the left in the picture, we can moreover deduce that the cells are steeper on the 
upflow side than on the downflow one. 
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FIG. 11. Deviation of the cell shape from a parabola, for a case without flow (solid line) and one with 
a flow / = 2.0 (dashed line). G = 0.5. 

V. CONCLUSIONS 

The addition of a shear flow to directional solidification near the absolute stabihty hniit affects 
the system in several more or less profound ways. 

First, the flow breaks the mirror symmetry of the equation of motion. This implies that cellular 
solutions are asymmetric in general. Symmetry breaking of this type is known to lead to drifting 
solutions in cases where it appears spontaneously |2|,Q. It produces the same behavior here. The 
appearance of a drift velocity can be understood at the level of a linear stability analysis which 
moreover provides a decent quantitative estimate for its value. As expected, this description of the 
drifting pattern breaks down for larger flows, where nonlinear effects exert a stronger influence. 

Second, the patterns appearing at the instability threshold of the planar interface for G > 
Gc{f = 0) are not hexagons but stripes. We shall see that this statement will have to be made 
more precise in II, because at this moment, we cannot say anything about the stability of stripe 
structures. The planar front becomes unstable via a supercritical Hopf bifurcation. Since the 
transition to hexagons of the system without flow is transcritical, i.e., hexagons can exist even 
below the threshold of the diffusive instability, we may expect stripes and hexagons to interact, at 
least at small flow strengths. This point will be discussed in some detail in II. 

The main effect of the flow in situations where a structure of hexagonal cells develops (i.e., for 
G < Gc{f = 0) and not too large a flow) is to promote ordering, i.e., the appearance of a hexagonal 
translational lattice. Hence, defects are eliminated more efficiently in the laterally moving pattern 
than in one at rest. Defects that stay tend to become aligned in the flow to give grain boundaries a 
preferential orientation (perpendicular to the flow). Finally, the flow increases instability towards 
local oscillations of cells with relative phase shifts. Since these oscillations also act to reduce the 
number of defects p^ , the flow reinforces this tendency, again working to improve translational 
order. 
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APPENDIX A: REAL SPACE REPRESENTATION OF THE FLOW TERM 

Since the nonlocal term is given by a product in Fourier space 

c[cr^\p\c, (Ai) 
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(where for brevity we denote the Fourier transform by an asterisk) its position space expression can 
be obtained as a convolution integral. Care must be taken, however, because the inverse Fourier 

transform of |p| docs not exist. Therefore, wc first rewrite €[(]* as a different product. 
For a one-dimensional interface, the most convenient approach seems to be to write 



jC[C]* = sign(p)pC* = -isign(p)(5xC)* 



(A2) 



While the (inverse) Fourier transform of the sign function does not exist as a function, it is defined 
in the distribution sense and easily calculable: 

s{x) = e^f-sign(p)dp = ^ e'f-rfp - e-'^^dp) = ^ | , (A3) 

the last expression being the distribution that is pointwise equal to i /ttx but requires any integral 
in which it appears to be interpreted as a principal value. 
Hence we obtain for the flow term in ID: 

f f°° 1 
-fd,C[C] = -^dA dx' ^d,<{x') , 

TT J _^ X-X' 

where the bar indicates a principal value integral. More specifically, we define 



(A4) 



dx' . . . = lim 

£-0+ 



dx' . 



dx' . 



For a two-dimensional interface, another decomposition is more appropriate 

c[cr 



^(-p^)c = -^(v^cr 



(A5) 



(A6) 



(In one dimension, the inverse Fourier transform of l/|p| is problematic because of the divergence 
at the origin, which in two dimensions is compensated by the volume element.) Hence we have 



/OO j*00 
dx' d2/'s(x-x')V'\(x'), 
-CSO J — OO 



(A7) 



where [p = {Px,Py)] 



1 r 



dpy 



-1 



47r2 



dp- exp (iplxl cos ^) . 
P 



(A8) 



To obtain the second expression, we have oriented the polar coordinate system such that the ray 
(/) = is parallel to x (hence px = p|x| cos^). Therefore, 



27r 



d0 



■ 7r(5(|x| cos ( 



47r2 |x| Jo cos (/) A-K Jq 



X sm( 



27r|x| 



(A9) 



The principal value integral of 1/ cos (f) vanishes as it extends over an entire period. Thus we arrive 
at the following final expression for the flow term 



-fo.m = i 9. y_ dx' y_ ^.'^ v'^c(x') . 



(AlO) 
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Both the ID and 2D expressions clearly exhibit the nonlocality of the flow term and its odd 
parity symmetry. All the other terms in Eq. (|lj) have even parity, so the flow term provides a 
symmetry breaking mechanism. 

Note also that the nonlocal kernels appearing in these expression are almost the simplest possible, 
if one thinks in terms of gradient expansions. Local terms produce successive powers of derivatives, 
corresponding to powers of p in Fourier space. Nonlocal terms would be connected with powers 
of l/|p|, which must be compensated for (in order to keep things finite at small p) by spatial 
derivatives. 



APPENDIX B: ANALYTIC CALCULATION OF DRIFT VELOCITIES 



In order to solve Eqs. and ( |2(]| ) analytically, we first cast them in a simpler form introducing 

UJr^UJr + 2q^ . (Bl) 

This yields 

Qt + CjI (8G -q"- 8g2) _ = o , (B2) 

fq^ {Cjr ~ 2q^) - 2(I;3 (-4^)^ ~ 2q^ - 8) = . (B3) 

We then consider the limiting cases / ^ 1 and / ^ 1. The first of these is very simple. As we 



know that there are solutions without flow, we set, as a first approximation, / = in Eqs. (B2) 



and (B3). uji will still remain / dependent via (iq), which takes the form 



UJi — 

2uJr 

We get immediately 

1 



(B4) 



Ur^2+ ^q^ (B5) 



and the Eq. (p3) for q can be reduced to quadratic. The result is 



1' 



(jJr 



-4+— V4 + 2G, (B6) 



-2q2+ V4 + 2G, (B7) 
v3 




c^, = / I 1 - ^== ) , (B8) 

(B9) 



1 2V4+2G 



As an example, for G — 0.7, this gives a drift velocity v — 0.218/, and the dependence on G is 
weak. Note that for this G value < and hence the structure will decay to a planar front 
(moving along as it does so at the the calculated velocity). 



The case / ^ 1 is slightly more complicated. Considering (B2) and the fact that G is of order 1 
in the interesting parameter range, we see (from the signs in the equation) that the only possible 
dominant balance for a pair of terms is 

^t--^- (BIO) 
Inserting this in (|B3|) and combining the two equations, we find however 
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(Bll) 



which is in contradiction with (BIO) for / 1 and q = 0(1). This imphcs that dominant balances 
in (B2) must contain three terms at least and that th e wa venum ber q of the fastest-growing mode 
has to scale, too. It must increase with / and Eqs. ( B10| ) and (Bll) suggest the scaling / ^ q^. 
Hence, we set 



f = aq^ , 

assuming a = 0{1). Inserting this into Eqs. (p32|) and (B3) and using 5 ^ 1, we arrive at 



~4 -2 4 a^?^ „ 
U}^ — ujj.q — = , 



Equation ( ^ ) can be solved for o)^: 



r,2 _ 1 



(1 + Vl + a') 



Using this in (B14), we end up, after some simplifications, with a relation determining a: 

~8" 



(1 + vi+a^)^^^ - yrr^ 



= . 



(B12) 

(B13) 
(B14) 

(B15) 
(B16) 



The numerical solution of this algebraic equation yields a — 28.88. (There is only one real solution.) 
We then obtain 



« 0.186 V7, 



1 + x/TT^ 



- 2 - « 0.065/ . 

a 



2/ 

-0.695 V7. 

q 



0.129/, 



(B17) 
(B18) 
(B19) 
(B20) 



All of these expressions are independent of G, as they must, being leading order results for / 3> G. 
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